单细胞转录组基础分析二：数据质控与标准化

单细胞测序技术是近年最大的生命科学突破之一，相关文章频繁发表于各大顶级期刊，然而单细胞数据的分析依然是大家普遍面临的障碍。本专题将针对10X Genomics单细胞转录组数据演示各种主流分析，包括基于Seurat的基础分析、以及基于clusterProfiler、Monocle、SingleR等R包的延伸分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨！（扫描文末二维码）

本节及之后的三、四节主要介绍单细胞表达矩阵到细胞类型鉴定的分析步骤，流程图如下：

10X官网下载数据

10X官网演示数据：

https://support.10xgenomics.com/single-cell-gene-expression/datasets

官方演示数据集的页面如下，Chromium Connect是自动化系统，官方介绍操作造成的批次效应较小。Chromium Next GEM最新版的芯片，由老版的双十字微流控芯片改成了单十字微流控芯片。目前国内的10X试剂基本都是V3版本的了，因此也不要下载V1和V2试剂的数据了。为了保证笔记本电脑能运行，我们下载细胞数比较少的数据，下图红色箭头所示：

本教程的分析都是基于箭头标示的数据。点击之后需要提交个人信息，提交之后就进入数据下载界面 了。

下载红框标示的数据，Seurat分析的输入文件在这里，解压后如下图所示：

数据质控与标准化 

上游分析软件Cell Ranger会对数据进行质控，但是在进行下游分析前，我们一般会对数据进行更严格的过滤。

library(Seurat)
library(tidyverse)
rm(list=ls())
dir.create("QC")
##创建seurat对象
scRNA.counts <- Read10X(data.dir = "filtered_feature_bc_matrix")    
try({scRNA = CreateSeuratObject(scRNA.counts[['Gene Expression']])},silent=T)
if(exists('scRNA')){} else {scRNA = CreateSeuratObject(scRNA.counts)}
#table([e-mail]$orig.ident)         #查看样本的细胞数量

##计算质控指标
#计算细胞中核糖体基因比例
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
#计算红细胞比例
HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB_m <- match(HB.genes, rownames(scRNA@assays$RNA)) 
HB.genes <- rownames(scRNA@assays$RNA)[HB_m] 
HB.genes <- HB.genes[!is.na(HB.genes)] 
scRNA[["percent.HB"]]<-PercentageFeatureSet(scRNA, features=HB.genes) 
#head([e-mail])
col.num <- length(levels([e-mail]))
violin <- VlnPlot(scRNA,
          features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"), 
          cols =rainbow(col.num), 
          pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
          ncol = 4) + 
          theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) 
ggsave("QC/vlnplot_before_qc.pdf", plot = violin, width = 12, height = 6) 
ggsave("QC/vlnplot_before_qc.png", plot = violin, width = 12, height = 6)  
plot1 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot3 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.HB")
pearplot <- CombinePlots(plots = list(plot1, plot2, plot3), nrow=1, legend="none") 
ggsave("QC/pearplot_before_qc.pdf", plot = pearplot, width = 12, height = 5) 
ggsave("QC/pearplot_before_qc.png", plot = pearplot, width = 12, height = 5)

运行上面的代码后会在"QC"文件夹下面得到4个文件

打vlnplot_before_qc的文件

上面的小提琴图依次是细胞的基因数量、mRNA分子数量、线粒体基因比例、红细胞基因比例。我们一般根据基因数量和线粒体比例来过滤细胞，细胞的最低基因数一般选择200-500，最大基因数和核糖体比例根据上图来选择，我的建议是minGene=500，maxGene=4000，pctMT=15。

##设置质控标准
print(c("请输入允许基因数和核糖体比例，示例如下：", "minGene=500", "maxGene=4000", "pctMT=20"))
minGene=500
maxGene=4000
pctMT=15

##数据质控
scRNA <- subset(scRNA, subset = nFeature_RNA > minGene & nFeature_RNA < maxGene & percent.mt < pctMT)
col.num <- length(levels([e-mail]))
violin <-VlnPlot(scRNA,
         features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"), 
         cols =rainbow(col.num), 
         pt.size = 0.1, 
         ncol = 4) + 
         theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) 
ggsave("QC/vlnplot_after_qc.pdf", plot = violin, width = 12, height = 6) 
ggsave("QC/vlnplot_after_qc.png", plot = violin, width = 12, height = 6)
scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize", scale.factor = 10000)

##保存中间结果
saveRDS(scRNA, file="scRNA.rds")

运行上面的代码后会在"QC"文件夹下面得到vlnplot_after_qc的文件

可以看到基因数和核糖体比例不正常的细胞都过滤了。以上代码的最后一行是对数据进行标准化，它的作用是让测序数据量不同的细胞的基因表达量具有可比性。计算公式如下：

标准化后基因表达量 = log1p（10000*基因counts/细胞总counts）

获取帮助

本教程的目的在于把常用的单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以分享此篇文章到朋友圈，截图微信发给Kinesin（文末二维码），我会抽时间组群直播讲解单细胞数据分析的全过程。本专题所用的软件、数据、代码脚本和中间结果，我打包放在了百度云上，需要的朋友添加Kinesin微信索取。